General approach to the assignment

Through this assignment, I want to see the 5-mins walkshed and the bikeshed of the polling locations in South Boston. I will create 5-mins isochrones for each transportation mode, and I will calculate and compare them. I will also calculate the number of trees in each walkshed of the polling locations in South Boston and visualize it. I’m a strong believer that fresh air makes a significant impact on a person’s smart decision!

library(dplyr)
library(knitr)
library(osmdata)
library(opentripplanner)
library(tidyverse)
library(sf)
library(ggthemes)
library(ggspatial)
library(splitstackshape)
MA_state_plane <- "+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs"

Downloading polling location data in Boston

poll_location <- st_read(
  "http://bostonopendata-boston.opendata.arcgis.com/datasets/f7c6dc9eb6b14463a3dd87451beba13f_5.kml?outSR=%7B%22latestWkid%22%3A2249%2C%22wkid%22%3A102686%7D")

Downloading Boston neighborhood data

boston_nhood <- st_read("http://bostonopendata-boston.opendata.arcgis.com/datasets/3525b0ee6e6b427f9aab5d0a1d0a1a28_0.kml?outSR=%7B%22latestWkid%22%3A2249%2C%22wkid%22%3A102686%7D") %>%
  filter(Name == "South Boston")

Filtering polling locations in South Boston

south_boston_polls <- poll_location[boston_nhood, ]

Get street data - South Boston

opq(bbox = 'South Boston MA USA') %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_xml(file = 'OTP/graphs/default/south_boston_street.osm')
south_boston_street_features <- opq(bbox = 'South Boston MA USA') %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_sf()

south_boston_streets <- south_boston_street_features$osm_lines %>%
  st_transform(crs = MA_state_plane)
ggplot(south_boston_streets) +
  geom_sf() +
  theme_map()

Setting up Open Trip Planner

path_data <- file.path(getwd(), "OTP")
path_otp <- paste(path_data, "otp.jar",sep = "/")

otp_build_graph(otp = path_otp, dir = path_data, memory = 1024) 
otp_setup(otp = path_otp, dir = path_data, memory =1024)
## 2020-10-07 12:41:54 OTP is loading and may take a while to be useable
## Router http://localhost:8080/otp/routers/default exists
## 2020-10-07 12:42:55 OTP is ready to use Go to localhost:8080 in your browser to view the OTP
otpcon <- otp_connect()
## Router http://localhost:8080/otp/routers/default exists

Filtering polling locations in South Boston

south_boston_polls <- poll_location[boston_nhood, ] %>%
  st_transform(src = MA_state_plane)

Creating isochrones

I’m going to create 5-mins isochrones for walking and bicycling.

iso_5min_walk <- 
  otp_isochrone(otpcon = otpcon, fromPlace = south_boston_polls, 
                mode = "WALK", cutoffSec = 300) %>%
  st_transform(crs = MA_state_plane) %>%
  mutate(mode = "walk")
## Warning in otp_isochrone(otpcon = otpcon, fromPlace = south_boston_polls, :
## Failed to get isochrone with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 300},"id":"fid-c895b39_174fe96b228_-7f9d"}]}Failed to get
## isochrone with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 300},"id":"fid-c895b39_174fe96b228_-7f96"}]}Failed to get
## isochrone with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":300},"id":"fid-
## c895b39_174fe96b228_-7f94"}]}
iso_5min_bike <- 
  otp_isochrone(otpcon = otpcon, fromPlace = south_boston_polls, 
                mode = "BICYCLE", cutoffSec = 300) %>%
  st_transform(crs = MA_state_plane) %>%
  mutate(mode = "bike")

iso_all_modes <- rbind(iso_5min_bike, iso_5min_walk)

Isochrones Map (1)

right_side <- st_bbox(iso_all_modes)$xmax
left_side  <- st_bbox(iso_all_modes)$xmin
top_side <- st_bbox(iso_all_modes)$ymax
bottom_side <- st_bbox(iso_all_modes)$ymin

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = south_boston_polls) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("By bike", "By foot")) +
  theme_map() +
  theme(legend.background = element_rect(fill = alpha("white", 0.5),
                                         color = "gray")) + 
  labs(caption = "Basemap Copyright OpenStreetMap contributors")
## Loading required namespace: raster

Isochrones Map (2)

right_side <- st_bbox(iso_all_modes)$xmax
left_side  <- st_bbox(iso_all_modes)$xmin
top_side <- st_bbox(iso_all_modes)$ymax
bottom_side <- st_bbox(iso_all_modes)$ymin

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type = "stamenbw", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = south_boston_polls) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("By bike", "By foot")) +
  theme_map() +
  theme(legend.background = element_rect(fill = alpha("white", 0.5),
                                         color = "gray")) + 
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

Calcuate and compare isochrone areas

south_boston_polls <- st_difference(south_boston_polls)
## although coordinates are longitude/latitude, st_difference assumes that they are planar
iso_areas <- iso_all_modes %>%
  mutate(area = st_area(iso_all_modes)) %>%
  st_set_geometry(NULL) %>%
  pivot_wider(names_from = mode, values_from = area) 
## Warning: Values are not uniquely identified; output will contain list-cols.
## * Use `values_fn = list` to suppress this warning.
## * Use `values_fn = length` to identify where the duplicates arise
## * Use `values_fn = {summary_fun}` to summarise duplicates
iso_areas <- iso_areas %>%
  filter(walk != "NULL") %>%
  filter(str_detect(walk,"c")==FALSE) %>%
  filter(str_detect(bike,"c")==FALSE)
## Warning in stri_detect_regex(string, pattern, negate = negate, opts_regex =
## opts(pattern)): argument is not an atomic vector; coercing
ggplot(iso_areas, 
       aes(x = as.numeric(walk), y = as.numeric(bike))) +
  geom_point() +
  scale_x_continuous(name = 
            "Area within a five-minute walking distance\nof a polling location(square km)",
            breaks = breaks <- seq(10000, 130000, by = 20000),
            labels = breaks / 1000000) +
  
  scale_y_continuous(name = 
            "Area within a five-minute biking distance\nof a polling location\n(square km)",
            breaks = breaks <- seq(0, 700000, by = 100000),
            labels = breaks / 1000000) +
  theme_bw()

Counting the number of trees within 5-min walk from polling locations in South Boston

tree <- st_read("http://bostonopendata-boston.opendata.arcgis.com/datasets/ce863d38db284efe83555caf8a832e2a_1.kml", quiet = TRUE) 
south_boston_tree <- tree[boston_nhood, ] 
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
ggplot(boston_nhood) +
  geom_sf() + 
  geom_sf(data = south_boston_tree)

south_boston_tree <- south_boston_tree %>%
  st_transform(MA_state_plane)
tree_poll_walk <- south_boston_tree [iso_5min_walk, ]

ggplot(tree_poll_walk) +
  geom_sf(data = iso_5min_walk, fill = "royalblue3", alpha = 0.5, color = NA) + 
  geom_sf(data = tree_poll_walk , color = "darkgreen",size = 0.3) +
  theme_map()

south_boston_tree <- south_boston_tree %>%
  st_join(tree_poll_walk) %>%
  mutate(by_tree = !is.na(Name.y))
n_tree_by_walk <- sum(south_boston_tree$by_tree)

n_tree_by_walk
## [1] 411
n_south_boston_tree <- length(south_boston_tree$by_tree)
pct_tree_by_walk <- n_tree_by_walk/n_south_boston_tree
pct_tree_by_walk
## [1] 0.06210335
left_side <- st_bbox(south_boston_tree)$xmin
top_side <- st_bbox(south_boston_tree)$ymax

ggplot(south_boston_tree) +
  geom_sf(fill = "lightgray", color = "darkgray") +
  geom_sf(data = south_boston_tree, size = 0.3,
          aes(color = by_tree)) +
  scale_color_manual(values = c("firebrick3", "darkgreen"),
                    name = "South Boston Trees \nby a distance to poll locations",
                    labels = c("No Tree within 5-min walk",
                               "Tree within 5-min walk")) +
  annotation_scale(location = "br") +
  annotation_north_arrow(location = "tr", 
                         style = north_arrow_minimal()) +
  annotate(geom = "text", x = right_side, y = bottom_side+200,
           label = paste("Of the ", 
                      prettyNum(n_south_boston_tree, big.mark = ","),
                       " trees in South Boston,\n", 
                       prettyNum(n_tree_by_walk, big.mark = ","),
                       " (", 
                       prettyNum(100*pct_tree_by_walk, digits = 0),
                       "%) are within \n5-mins walk from polling locations",
                       sep = ""), size = 3) +
  theme_map() +
  theme(panel.background = element_rect(fill = "aliceblue"),
        legend.background = element_rect(fill = alpha("white", 0.5), 
                                         color = "gray"))